package scu.maqiang.homogenization;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L21;
import scu.maqiang.fes.FES1L2PE;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class PESOTS1D {

	public static void main(String[] args) {
		long begin = System.nanoTime();
		int N = 32;
		double ep = 1.0 / N;
		Mesh1L2 cellMesh = new Mesh1L2();
		cellMesh.line(0.0, 1.0, 128);
		ScalarFunc regionCellFunc1D = (x, llable, param) -> {
			return Math.abs(x[0] - 0.5) < 0.25 ? 1.0 : -1.0;
		};
		cellMesh.setDomainLabel(regionCellFunc1D, null, 10);
		FES1L2PE fsCell = new FES1L2PE(cellMesh);
		double[][] FOCS = new double[2][fsCell.GetNdof()];
		double[][] SOCS = new double[3][fsCell.GetNdof()];
		double[][] D = new double[2][3];
		SOTS1DL2PE ss = new SOTS1DL2PE(fsCell);
		ScalarFunc CellG1 = (xyz, label, param) -> label == 10 ? 20 : 4;
		ScalarFunc CellEp1 = (xyz, label, param) -> label == 10 ? 1.0 : 0.2;
		ScalarFunc CellE1 = (xyz, label, param) -> label == 10 ? 5.0 : 3.0;
		ScalarFunc Cellrho = (xyz, label, param) -> label == 10 ? 1 : 1;
		ScalarFunc Cellf1 = (xyz, label, param) -> label == 10 ? 1 : 1;
		ScalarFunc Cellrf = (xyz, label, param) -> label == 10 ? 1 : 1;
		ss.PECellSolution(new ScalarFunc[] { CellG1, CellEp1, CellE1, Cellrho, Cellf1, Cellrf }, null, FOCS, D,
				SOCS);
		double[][] FOCSUV1 = new double[2][cellMesh.getNv()];
		double[][] FOCSUV2 = new double[2][cellMesh.getNv()];
		fsCell.extractComponent(FOCS[0], FOCSUV1);
		fsCell.extractComponent(FOCS[1], FOCSUV2);
		double[][] FFOCS = new double[4][];
		FFOCS[0] = FOCSUV1[0];
		FFOCS[1] = FOCSUV1[1];
		FFOCS[2] = FOCSUV2[0];
		FFOCS[3] = FOCSUV2[1];
		System.out.println("N1 min: " + MVO.min(FOCSUV1[0]) + "\tmax: " + MVO.max(FOCSUV1[0]));
		System.out.println("H  min: " + MVO.min(FOCSUV1[1]) + "\tmax: " + MVO.max(FOCSUV1[1]));
		System.out.println("F  min: " + MVO.min(FOCSUV2[0]) + "\tmax: " + MVO.max(FOCSUV2[0]));
		System.out.println("R1 min: " + MVO.min(FOCSUV2[1]) + "\tmax: " + MVO.max(FOCSUV2[1]));

		cellMesh.toTecplot("CellFOCS.dat", FFOCS);
		// System.out.println("min: " + MVO.min(FOCS[0]) + "\tmax: " +
		// MVO.max(FOCS[0]));
		// System.out.println("min: " + MVO.min(FOCS[1]) + "\tmax: " +
		// MVO.max(FOCS[1]));
		System.out.println("Homogenized PE Coefficients: ");
		System.out.println(MVO.toString(D));
		ss.computeHomogenizedPE(new double[] { 20, 1.0, 5.0 }, new double[] { 4, 0.2, 3.0 }, 0.5);
		System.out.println("min: " + MVO.min(SOCS[0]) + "\tmax: " + MVO.max(SOCS[0]));

		double[][] SOCSUV1 = new double[2][cellMesh.getNv()];
		double[][] SOCSUV2 = new double[2][cellMesh.getNv()];
		double[][] SOCSUV3 = new double[2][cellMesh.getNv()];
		fsCell.extractComponent(SOCS[0], SOCSUV1);
		fsCell.extractComponent(SOCS[1], SOCSUV2);
		fsCell.extractComponent(SOCS[2], SOCSUV3);
		double[][] SSOCS = new double[6][];
		SSOCS[0] = SOCSUV1[0];
		SSOCS[1] = SOCSUV1[1];
		SSOCS[2] = SOCSUV2[0];
		SSOCS[3] = SOCSUV2[1];
		SSOCS[4] = SOCSUV3[0];
		SSOCS[5] = SOCSUV3[1];
		cellMesh.toTecplot("CellSOCS.dat", SSOCS);

		// ���Ȼ����
		Mesh1L2 homoMesh = new Mesh1L2();
		homoMesh.line(0, 1, 256);
		FES1L2PE fsHomo = new FES1L2PE(homoMesh);
		SRMatrix homoA = new SRMatrix(fsHomo.GetNdof());
		double[] homoSol = new double[fsHomo.GetNdof()];
		double[] homoRHS = new double[fsHomo.GetNdof()];
		fsHomo.assemblePEStiff(new double[] { D[0][0], D[0][1], D[0][2] }, BVPType.COMMON, homoA);
		fsHomo.assemblePESource(new double[] { D[1][1], D[1][2] }, BVPType.COMMON, homoRHS);
		fsHomo.applyBC_MBN(homoA, Direct.All, 1);
		fsHomo.applyBC_MBN(homoRHS, Direct.All, 0.0, 1);
		DirectSSolver homoSolver = new DirectSSolver(homoA);
		homoSol = homoSolver.LUSolve(homoRHS, 1, 1e-6);

		Mesh1L2 fineMesh = new Mesh1L2();
		fineMesh.line(0, 1, 32768);
		ScalarFunc regionFineFunc1D = (x, llable, param) -> {
			double xCell = x[0] * N - (int) (x[0] * N);
			return Math.abs(xCell - 0.5) < 0.25 ? 1.0 : -1.0;
		};
		fineMesh.setDomainLabel(regionFineFunc1D, null, 10);
		FES1L2PE fsFine = new FES1L2PE(fineMesh);
		SRMatrix fineA = new SRMatrix(fsFine.GetNdof());
		double[] fineSol = new double[fsFine.GetNdof()];
		double[] fineRHS = new double[fsFine.GetNdof()];
		fsFine.assemblePEStiff(new ScalarFunc[] { CellG1, CellEp1, CellE1 }, null, BVPType.COMMON, fineA);
		fsFine.assemblePESource(new ScalarFunc[] { Cellrho, Cellrf }, null, BVPType.COMMON, fineRHS);
		fsFine.applyBC_MBN(fineA, Direct.All, 1);
		fsFine.applyBC_MBN(fineRHS, Direct.All, 0.0, 1);
		DirectSSolver fineSolver = new DirectSSolver(fineA);
		fineSol = fineSolver.LUSolve(fineRHS, 1, 0.005);

		double[] HomoFineSol = new double[fsFine.GetNdof()];
		double[] FirstSol = new double[fsFine.GetNdof()];
		double[] SecondSol = new double[fsFine.GetNdof()];
		ss.assemblePESOTSSolution(fsFine, HomoFineSol, FirstSol, SecondSol, fsHomo, homoSol, fsCell, FOCS, SOCS, N, ep);

		System.out.println("Homo Solution");
		System.out.println("Homo min: " + MVO.min(HomoFineSol) + "\tmax: " + MVO.max(HomoFineSol));
		System.out.println("First Order Approximation: ");
		System.out.println("FOA min: " + MVO.min(FirstSol) + "\tmax: " + MVO.max(FirstSol));
		System.out.println("Second Order Approximation: ");
		System.out.println("SOA min: " + MVO.min(SecondSol) + "\tmax: " + MVO.max(SecondSol));

		double[][] AllSol = new double[8][fineMesh.getNv()];
		double[][] HomoSolUV = new double[2][fineMesh.getNv()];
		fsFine.extractComponent(HomoFineSol, HomoSolUV);
		double[][] FirstSolUV = new double[2][fineMesh.getNv()];
		fsFine.extractComponent(FirstSol, FirstSolUV);
		double[][] SecondSolUV = new double[2][fineMesh.getNv()];
		fsFine.extractComponent(SecondSol, SecondSolUV);
		double[][] FineSolUV = new double[2][fineMesh.getNv()];
		fsFine.extractComponent(fineSol, FineSolUV);

		AllSol[0] = FineSolUV[0];
		AllSol[1] = FineSolUV[1];
		AllSol[2] = HomoSolUV[0];
		AllSol[3] = HomoSolUV[1];
		AllSol[4] = FirstSolUV[0];
		AllSol[5] = FirstSolUV[1];
		AllSol[6] = SecondSolUV[0];
		AllSol[7] = SecondSolUV[1];
		fineMesh.toTecplot("PESOTSSolution.dat", AllSol);

		FES1L21 fsL2 = new FES1L21(fineMesh);

		double L2Norm = fsL2.computeL2Norm(FineSolUV[0]);
		double[] error0 = MVO.add(FineSolUV[0], -1.0, HomoSolUV[0]);
		double[] error1 = MVO.add(FineSolUV[0], -1.0, FirstSolUV[0]);
		double[] error2 = MVO.add(FineSolUV[0], -1.0, SecondSolUV[0]);
		System.out.println("Displacement L2Rrror:");
		System.out.println("Error 0: " + fsL2.computeL2Norm(error0) / L2Norm);
		System.out.println("Error 1: " + fsL2.computeL2Norm(error1) / L2Norm);
		System.out.println("Error 2: " + fsL2.computeL2Norm(error2) / L2Norm);

		double H1SemiNorm = fsL2.computeH1SemiNorm(FineSolUV[0]);
		System.out.println("Displacement H1 SemiRrror:");
		System.out.println("Error 0: " + fsL2.computeH1SemiNorm(error0) / H1SemiNorm);
		System.out.println("Error 1: " + fsL2.computeH1SemiNorm(error1) / H1SemiNorm);
		System.out.println("Error 2: " + fsL2.computeH1SemiNorm(error2) / H1SemiNorm);
		System.out.println();
		
		L2Norm = fsL2.computeL2Norm(FineSolUV[1]);
		error0 = MVO.add(FineSolUV[1], -1.0, HomoSolUV[1]);
		error1 = MVO.add(FineSolUV[1], -1.0, FirstSolUV[1]);
		error2 = MVO.add(FineSolUV[1], -1.0, SecondSolUV[1]);
		System.out.println("Electric potential L2Rrror:");
		System.out.println("Error 0: " + fsL2.computeL2Norm(error0) / L2Norm);
		System.out.println("Error 1: " + fsL2.computeL2Norm(error1) / L2Norm);
		System.out.println("Error 2: " + fsL2.computeL2Norm(error2) / L2Norm);

		H1SemiNorm = fsL2.computeH1SemiNorm(FineSolUV[1]);
		System.out.println("Electric potential H1 SemiRrror:");
		System.out.println("Error 0: " + fsL2.computeH1SemiNorm(error0) / H1SemiNorm);
		System.out.println("Error 1: " + fsL2.computeH1SemiNorm(error1) / H1SemiNorm);
		System.out.println("Error 2: " + fsL2.computeH1SemiNorm(error2) / H1SemiNorm);
		long end = System.nanoTime();
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
	}

}
